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Abstract. We explore systematic biases in the identification of dark matter in future direct 
detection experiments and compare the reconstructed dark matter properties when assum- 
ing a self-consistent dark matter distribution function and the standard Maxwellian velocity 
distribution. We find that the systematic bias on the dark matter mass and cross-section de- 
termination arising from wrong assumptions for its distribution function is of order ~ la. A 
much larger systematic bias can arise if wrong assumptions are made on the underlying Milky 
Way mass model. However, in both cases the bias is substantially mitigated by marginalizing 
over galactic model parameters. We additionally show that the velocity distribution can be 
reconstructed in an unbiased manner for typical dark matter parameters. Our results high- 
light both the robustness of the dark matter mass and cross-section determination using the 
standard Maxwellian velocity distribution and the importance of accounting for astrophysical 
uncertainties in a statistically consistent fashion. 
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1 Introduction 

Determining the properties of weakly interacting massive particles (WIMPs) from the va- 
riety of available experimental techniques presents an exciting experimental and theoretical 
challenge [1-5]. Direct, indirect, and collider dark matter searches each come with their own 
unique set of benefits from a theoretical perspective, as well as systematic issues from an ex- 
perimental perspective. A solution to the dark matter problem certainly requires convergence 
from each of these major techniques (see e.g. Refs. [6-14]). 

This paper focuses on understanding how accurately WIMP properties - in particular 
the mass and elastic spin-independent (SI) scattering cross-section - can be determined from 
direct dark matter detection experiments given systematic uncertainties coming from our 
imperfect knowledge of the dark matter astrophysical distribution. Several recent studies have 
tackled this issue in one form or another: Strigari & Trotta [15] introduced a Bayesian method 
to determine WIMP properties by using kinematic data in the Milky Way and marginalizing 
over galactic halo model parameters; this method was then improved in Ref. [10] to consider 
complementarity between different experimental targets for a more general parameterization 
of the velocity distribution [16]. Using mock data from future direct detection experiments, 
Peter [17] (see also [18]) identified biases in the WIMP mass determination for different 
phenomenological velocity distribution models, determining both the circular and escape 
velocities from the event rate distribution, while [19] investigated systematic biases in the 
determination of WIMP properties arising from unavoidable statistical fluctuations in the 
energy distribution. 

Here we extend the above studies and consider the bias in the reconstructed WIMP 
properties, improving on previous work in two specific ways. First, we study bias for several 
theoretically well-motivated models of the smooth WIMP velocity distribution. In particular, 
we apply the above techniques to theoretical models that establish a direct correspondence 
between the velocity distribution and the dark matter density profile in a self-consistent 
manner. While there have been several works that have explored phenomenologically how 
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these models impact the observed event rate [16, 20-23], we examine for the first time the 
bias in the reconstruction of WIMP properties for these models (see also [24]). Second, 
we study the bias in the WIMP mass and cross-section determination while simultaneously 
fitting kinematic data that probe the dark matter distribution of the Galaxy. This approach 
accounts for the effect of uncertainties in galactic model parameters in a statistically robust 
fashion, and we argue that this or a similar technique will be necessary in the case of a 
confirmed WIMP signal. 

We envision that the approach undertaken in this paper extends beyond the "discovery" 
phase of dark matter detection, anticipating a period when robust statistical procedures will 
be required to extract WIMP properties in an accurate way from direct detection data. Our 
approach is complementary to recent ones [25-27] that have been developed to elucidate 
whether the events seen by multiple experiments [28-31] can be consistent with a WIMP 
signal. These approaches encapsulate the impact of the astrophysical dark matter distribution 
in one single galactic halo model parameter that determines the event rate distribution and 
are of particular interest if several experiments are indeed seeing WIMP signals. 

The main results of this work quantify the bias that is incurred for different models of 
the velocity distribution, both with and without including uncertainties in galactic model pa- 
rameters. While we show that the WIMP mass and cross-section are systematically unbiased 
when assuming the true underlying velocity distribution (as one would expect), standard 
reconstruction methods do not recover the true event rate distribution when assuming an 
incorrect model of the velocity distribution. Nevertheless, we find that it is straightforward 
to overcome this bias by marginalizing over galactic model parameters, and that the use of 
the standard Maxwellian velocity distribution leads to a robust determination of the WIMP 
properties from direct detection experiments. 

The outline of the paper is as follows. In Section 2 we present the modeling of the 
different Milky Way components, while Section 3 details the calculation of the WIMP event 
rates. Section 4 is devoted to outlining our methodology. Our main results are shown and 
discussed in Section 5, before concluding in Section 6. 

2 Milky Way modelling 

In this section we motivate the models for the Milky Way that we consider, and relate these 
to the corresponding WIMP velocity distributions. 

2.1 Standard halo model 

For our main analysis we consider spherically symmetric and isotropic dynamical models. 
Though it is straightforward to consider more complicated models, the above assumptions 
are appropriate for the interpretation of the primary results of this paper. Let us define v 
as the WIMP velocity in the Earth rest frame, as the Earth velocity in the galactic rest 
frame, and w = v + vq as the WIMP velocity in the galactic rest frame. The Maxwellian 
distribution truncated at the escape velocity is 



for w < v csc and otherwise, with given values for and v esc . This velocity distribution, 
or its variations, are often referred to as the standard halo model (SHM). We take v§ ~ vq 




(2.1) 
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and define vo as the local circular velocity, which is derived from the total galactic potential, 
(hat, as 

(2.2) 



2 „ "Wot 
V ° = ~dr~ 

ar r=R e 

where Rq is the galactocentric distance. Similarly, the escape velocity is derived from the 
total potential as 

v 2 esc = 2\hot(R Q )\ . (2.3) 

To determine hot we consider a Milky Way model split into three separate components, the 
bulge (b), the disk (d) and the dark matter halo (dm), and sum the contribution from each 
component, 

hot =<t>b + <t>d + 4>dm ■ (2.4) 

For the bulge we take the spherical potential h = —GMi/(r + cq) and for the disk we take 
the spherical approximation hi = — GMd[l — exp(— r/bd)]/r, where Mb = 1.5 x 10 10 Mq is 
the fixed bulge mass, Co = 0.6 kpc is the fixed bulge scale radius, is the disk mass and bd 
is the disk scale radius (both treated as a free parameters, as explained later), see Strigari & 
Trotta [15] for references and further details. We specify the dark matter component by its 
mass density that is taken to follow a generalized Navarro-Frenk- White (NFW) profile: 

Pdm = (r/r s )«(l + (r/r s ) b )( c - a )/ b ' (2 ' 5) 

where a, b and c are profile indices, p s (r s ) is the scale density (radius) and po = pdm(R&) is 
the local dark matter density. 

With the potential specified, the radial velocity dispersion for a given stellar population 
is determined from the spherical and isotropic Jeans equation, 

%^ = - Ga ^ . (2.6) 
or dr 

where p* is the density distribution for the stellar population. As motivated in more detail 
below, we will assume that the population of stars is isotropic in velocity space. In the 
limit of an isotropic velocity distribution, the radial velocity dispersion a r from Eq. (2.6) is 
identical to the observed line-of-sight velocity dispersion. 

In this model, in order to evaluate Eq. (2.1), there are thus two parameters that must be 
determined from the kinematic data on the Galaxy, namely an d v esc . These are functions of 
the entire set of parameters that determine the total local potential, i.e. {p s , r s , a, b, c, bj, Mj, Rq}. 
These parameters are fit to the kinematic data described below. Once this is done, the WIMP 
event rate distribution is derived from the velocity distribution function in Eq. (2.1). 

The modeling described above is very similar to that used in Strigari & Trotta [15], 
who connected the WIMP event rate distribution to the kinematics of halo stars and other 
measurements of galactic model parameters. The only difference between what we have 
outlined above and the analysis by Strigari <fe Trotta is that here we set the anisotropy 
parameter (5 to zero. However, in this paper we extend the former formalism to study a self- 
consistent model that connects the dark matter density distribution with the corresponding 
velocity distribution. We now proceed to describe this model. 
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Figure 1. The WIMP velocity distribution (left) and the mean inverse velocity (right) for the case of 
the standard halo model (dashed lines) and the self-consistent distribution function (solid lines) with 
the fiducial values in Table 1. The mean inverse velocity -Fkstro m the right frame corresponds to the 
velocity integral in Eq. (3.1) for a 50 GeV WIMP and a xenon target. 



2.2 Self-consistent distribution function 

Though the above procedure for determining the velocity distribution from vq and v esc is 
well denned, it is not necessarily dynamically self-consistent, in that the velocity distribution 
in Eq. (2.1) does not follow from the bulge, disk, and dark matter potentials. It is thus 
important to compare the SHM to a self-consistent model, in which the distribution function 
is generated from an one-to-one mapping to the dark matter density profile. To perform this 
comparison, we consider the Eddington formula 

f (e) - — f &2pdm d * tot (2 7) 

which is valid under the assumption of spherical symmetry and isotropy. Here, e is the 
relative specific energy of the WIMP and reads 

e = -E(Rq) + &ot(oo) = -w 2 /2 + tttot(flo) (2-8) 

with E(Rq) = w 2 /2 + 0tot(-R©) and * to t(-Ro) = -</>tot(-R©) + ^tot(oo). For an assumed 
parameter set that determines the potential, the self-consistent distribution function (SCDF) 
derived from Eq. (2.7) differs from the SHM distribution. The crux of our work lies precisely 
in the comparison between the velocity distributions obtained from Eqs. (2.7) and (2.1) - 
both are plotted in Fig. 1 (left). 

Though deriving the velocity distribution self-consistently from the density profile is in 
principle more theoretically appealing than simply using the SHM, a few caveats are in order. 
First, it appears that the velocity distribution from cosmological simulations of dark matter 
halos differs significantly from the correspondence from Eddington theory [32-34]. Hence, 
this model itself may still be an oversimplification if the phase space properties of simulated 
halos are a reflection of the phase space properties of the Milky Way. Nonetheless, the 
comparison between the SHM and the SCDF has important phenomenological implications, 
as we show in detail below. 
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Table 1. Summary of the parameters used in the analysis. We have assumed uniform priors within 
the given ranges. The letter in the last column refers to the setups in which these parameters are 
used: f (fixed), p (pheno), and F (full). The setups are defined in Section 4. 



Symbol 


Parameter 


Fiducial values 


Prior 


range 


Model used 


logio i™x/GeVJ 


WIMP mass 


log 10 (25,50,250) 


0.1 - 


-3.0 


f,P,F 


logio kfVpbl 


spin-independent cross-section 


-9,-8.5 


-11 - 


-(-7) 


f,pF 


po |GeV/cm J j 


local dark matter density 


0.29 


0.001 


-0.79 


p 


Vq [km/sl 


circular velocity 


216 


66- 


366 


p 


^csc 


[km/s] 


escape velocity 


553 


388 - 


- 718 


p 


log 10 [p,/(M /kpc 3 )J 


scale density 


6.8 


5- 


-8 


F 


r s [kpc] 


scale radius 


20 


1 - 


60 


F 


a 




inner slope 


1 


0- 


- 2 


F 


b 




transition slope 


1 


0- 


- 2 


F 


c 




outer slope 


3 


2 - 


-5 


F 


b d [kpc] 


disk scale radius 


4 


2 - 


-6 


F 


M d 


10 10 M ] 


disk mass 


7 


2 - 


-8 


F 


Rq 


kpc] 


galactocentric distance 


8.33 


7- 


10 


P,F 



3 Direct Detection 



Our reconstruction of dark matter properties from simulated direct detection data follows 
standard analysis procedures; hence we only briefly outline the relevant points in this Section. 
The crucial galactic ingredient for direct detection is the dark matter phase space distribution, 
more precisely the local dark matter density po and local velocity distribution f(w). In fact, 
the scattering rate induced by WIMPs of mass m x and SI WIMP-proton cross-section a p l 
on target nuclei N(A, Z) is given by 

dR p Vp l , 2p2 ,„ , f d 3 



dER 2m x //| 



A 2 Fl{E R ) J^lf(w)e(v-v min ) , (3.1) 



where p p is the WIMP-proton reduced mass, v m \ n = (tunEr/ l {2p? N )) 1 / 2 is the minimum 
velocity of the WIMP in the Earth frame to produce a nuclear recoil of energy En and F$i 
is the spin-independent form factor (see [10] for further details). Eq. (3.1) assumes elastic 
Sl-nucleon scattering only; although other types of scattering have been extensively discussed 
in the literature, we focus on the elastic SI case for concreteness. 



4 Methodology 

4.1 Likelihood and data 

We adopt a Bayesian methodology and employ the MultiNest code [35-37] to derive samples 
from the posterior probability distribution function (pdf) p(Q\d) for the parameter set G 
given the data d using Bayes' theorem: 

p(G\d) oc £(0)p(9) (4.1) 

with p(0) the prior pdf (which is taken to be uniform within a certain range) and £(@) 
the likelihood function. In the case at hand, includes both the astrophysical parameters 
discussed in Section 2 and the WIMP properties m x and a p l - see Table 1, where the adopted 
prior ranges and fiducial values for each parameter are shown. 

We consider three separate setups, which differ in the number of parameters that are 
included in 0, the data sets used, and the models for the velocity distribution used: 
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• The first and simplest setup, which we label the fixed model (f), varies only the WIMP 
mass and cross-section, hence 

e = {m x ,a p SI } . (4.2) 

The galactic model parameters are fixed to their given fiducial values, whose motivation is 
discussed in more detail below. 

• The second setup, which we refer to as the phenomenological (i.e. pheno, p) model, 
varies the local dark matter density, local circular velocity, escape velocity, and galactocentric 
distance in addition to the WIMP mass and cross-section, so that 

® = { m x> ^p 1 ' Po, Vo,V csc , Rq} . (4.3) 

In this setup, we use the measurements of the local dark matter density [38-40], circular ve- 
locity [10, 41-45], escape velocity [46], and galactocentric distance [47], and adopt a Gaussian 
likelihood with the following means and standard deviations: 

po = 0.29 ±0.1 GeV/cm 3 
v = 216 ± 30 km/s 
^esc = 544 ± 33 km/s 

Rq = 8.33 ± 0.35 kpc , (4.4) 

where the central values for po and vq are the ones obtained in our Milky Way fiducial model, 
see below and Table 1. 

• The third setup (full, F) allows the largest flexibility in the galactic model parameters, 
as defined in Section 2, by treating all of them as free parameters along with the WIMP mass 
and cross-section. In this setup, the parameter set is 

= {m x , dp 1 , p s , r s , a, b, c, b d , M d , R Q } . (4.5) 

In addition to the above constraints on v esc and Rq, we add to the likelihood the constraint on 
the velocity dispersion profile of the Milky Way halo, utilizing the radial velocity dispersion 
of 2400 blue horizontal branch (BHB) halo stars as measured by the Sloan Digital Sky Survey 
(SDSS) [41]. In the isotropic Jeans equation, we take the stellar density profile of the tracer 
population to scale as p+ oc r~ 3 5 for r > 10 kpc and p± oc const otherwise. To implement 
Eq. (2.6), we use the measured velocity dispersion by Xue et al [41] into 9 radial bins, and 
assume an independent Gaussian likelihood for each bin. 

Note that, in the full setup, po, Vo, and v esc are derived parameters, and do not appear 
as an input into either the model for the smooth density distribution or in the velocity 
distribution. This is in contrast to the pheno setup, in which these parameters directly 
appear. In both the full and the pheno setups, po, vq, and v esc are constrained from the 
observations. 

The joint likelihood function £j mt(@) is the product of a term associated with galactic 
parameters and a term coming from direct detection data: 

•£joint(@) = ^astro 

(6)x£ DD (G) . (4.6) 
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The likelihood for the galactic parameters, £ a stro) is a multinormal Gaussian as described 
above. The direct detection likelihood, £dd> is computed as follows. First, we determine 
the fiducial values of the astrophysical parameters by performing a scan of the full model, 
using for the astrophysics likelihood the constraints for R@, v esc as well as the kinematics of 
SDSS halo stars discussed above. The mean of the posterior distribution over the astrophys- 
ical parameters is then used to approximately fix the fiducial values for each astrophysical 
parameter given in Table 1. This defines our assumed true galactic model (i.e., the fiducial 
model), making sure that it is close to what current astrophysical data require. We then 
generate mock direct detection counts through Eqs. (3.1) and (2.7) from the fiducial model, 
assuming a xenon, germanium or neon experiment with 2.00, 2.16 and 150 tonxyr of effec- 
tive exposures respectively (see Refs. [10, 48]) and several benchmark WIMP parameters. 
The energy resolution for the xenon and germanium cases is as in [10], while for neon we 
take the energy-dependent resolution a{En) = 1.0 keV^ EfkeV. The counts are distributed 
into 10 linearly-spaced bins in the range En = 10 — 100 keV for xenon and germanium and 
En = 20 — 100 keV for neon, and £dd is a standard binned likelihood. Since we are inter- 
ested in studying systematic biases arising from assumptions made about the dark matter 
distribution function, we generate mock counts without Poisson noise (rounding mock counts 
to the nearest integer, an approximation that is unlikely to have an impact on our results 
given the large number of counts expected in most energy bins). This procedure eliminates 
statistical fluctuations from our realised data, so that - in absence of systematic errors - our 
maximum a posteriori estimate of the parameters value ought to coincide with the true input 
value. 

4.2 Systematic bias from astrophysical assumptions 

With the ingredients outlined in the previous paragraphs we aim to quantify the system- 
atic bias caused by wrong assumptions for the galactic halo model. We thus run several 
scans attempting a reconstruction of the dark matter mass and cross-section from the above- 
described direct detection data, which were simulated from an astrophysical model obeying 
the Eddington formula (2.7), but reconstructing the WIMP properties adopting the truncated 
Maxwellian distribution (2.1). If the Eddington formula (2.7) and the truncated Maxwellian 
distribution (2.1) yield similar counts, or if astrophysical and statistical uncertainties are 
large, the posterior distribution should be well-centered about the true WIMP properties 
and the corresponding systematic bias is negligible. However, if that is not the case, then 
making the wrong assumption about the velocity distribution, i.e. adopting Eq. (2.1), will 
systematically bias the inference from direct detection data. 

In order to quantify this effect we define our "bias factor" a as 

with di t 2 = {( x t, Vt) — (x, y))-v\fl/\vi^\, where v\ and V2 {v\ and t^) are the eigenvalues (eigen- 
vectors) of the covariance matrix in the parameter space (x, y) = (log 10 [m x /GeV] , log 10 [crj^/pb] ), 
(xt,yt) indicates the true WIMP properties and (x,y) the posterior mean. The bias factor 
a simply measures the distance between the true and mean points in units of standard 
deviations along the eigenvectors of the covariance matrix. 



(4.7) 
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Figure 2. Unbiased reconstruction of WIMP properties in a xenon ton-scale instrument for the fixed 
astrophysical setup, when the distribution function used for the reconstruction ("rec") matches the 
one used to generate the data ("true"). Inner and outer contours correspond to the joint 68% and 
95% posterior probability contours, respectively. The red (dark) contours assume the standard halo 
model both in generating the mock data and in the reconstruction, while the green (light) contours 
assume the self-consistent distribution function. 



5 Results & Discussion 

Firstly, we investigate the statistical accuracy with which the WIMP mass and cross-section 
can be determined in the absence of systematic bias, i.e., if we assume that we know the 
correct velocity distribution model. The reconstruction is shown in Fig. 2, for three bench- 
mark values of the WIMP mass. In each case, we simulate counts assuming either the SHM 
or the SCDF, and we reconstruct mass and cross-section assuming the correct distribution 
model and correct, fixed parameters for it. This is of course an overly optimistic scenario, 
but it delineates the maximum statistical accuracy that can be achieved in the absence of 
astrophysical uncertainties. As shown by Fig. 2, the determination of the WIMP mass and 
cross-section is roughly equivalent for both the SHM and SCDF in this optimistic scenario. 
For each fiducial set of parameters, there is no significant systematic bias, i.e. the fiducial 
parameters are well-determined at the one-sigma level. 

We now study the systematic bias arising from assuming the wrong velocity distribution 
model in the reconstruction - for reference, the integral in Eq. (3.1) (which controls the impact 
of the different velocity distribution models on the scattering rate) is shown in Fig. 1 (right) in 
the case of the SHM and the SCDF. Fig. 3 shows an example of the bias that can be incurred 
by assuming the wrong model for the velocity distribution. In this case, the underlying true 
event rate distribution is simulated from the SCDF, while the reconstruction is done using 
the SHM velocity distribution. The galactic parameters are fixed to their fiducial values in 
Table 1. Fig. 3 indicates a moderate bias in the reconstruction (typically around the la 
level, as quantified by the bias parameter a) for dark matter masses below 100 GeV. This 
bias reduces significantly for WIMP masses larger than 100 GeV, but this is a consequence 
of the statistical errors becoming much larger at those mass values, as the mass-cross-section 
degeneracy means that only a lower bound can be placed on the WIMP mass. 

Fig. 4 shows the reconstruction for three WIMP benchmarks, but expands on the anal- 
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Figure 3. Bias in the reconstruction of WIMP properties in a xenon ton-scale instrument for the 
fixed astrophysical setup, when the reconstruction assumes the standard halo model while the data 
are generated from the self-consistent distribution function with the fiducial parameters specified in 
Table 1. Inner and outer contours correspond to the joint 68% and 95% posterior probability contours, 
respectively. The parameter a quantifies the systematic bias according to Eq. (4.7). 
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Figure 4. Reconstruction of WIMP properties when (wrongly) assuming the standard halo model 
while the data are generated from the self-consistent distribution function, but additionally marginal- 
izing over astrophysical uncertainties. The left frame shows the case of a xenon ton-scale instrument 
for the pheno and full astrophysical setups, while the right frame corresponds to the combination of 
xenon, germanium and neon data for the full astrophysical setup. Inner and outer contours corre- 
spond to the joint 68% and 95% posterior probability contours, respectively. The systematic bias, as 
quantified by a, is reduced with respect to Fig. 3, but statistical uncertainties increase. 



ysis of Fig. 3 by further marginalizing over galactic parameters. This is a more realistic 
scenario, where uncertainties in the astrophysical parameters are fully included in the in- 
ference about the WIMP properties by marginalizing over them. Encouragingly, we find 
that there is no significant bias in the reconstruction in this case even when assuming the 
wrong velocity distribution model. Thus, for the models of the velocity distribution that 
we compare, assuming the wrong velocity distribution does not significantly bias the mass 
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Figure 5. Reconstruction of WIMP properties in a xenon ton-scale instrument when assuming a 
wrong outer slope c and scale radius r s . Red (dark) contours are for data generated using the self- 
consistent distribution function with c = 3 and r s — 20 kpc, while the reconstruction wrongly assumes 
c = 4 and r s = 30 kpc (all other galactic parameters are kept fixed). Green (light) contours result 
from the same procedure but additionally marginalizing over the parameters c € [2, 5] and r s € [1, 60] 
kpc (with flat priors), thus significantly reducing systematic bias, as quantified by the parameter a. 

reconstruction. Of course, the price to pay for reducing the systematic bias is an increased 
statistical error on the WIMP parameters (compare with Fig. 3). Furthermore, we notice 
that marginalizing over astrophysical uncertainties in the more restricted pheno setup is suffi- 
cient to largely eliminate systematic bias, although the corresponding statistical uncertainty 
is higher than in the full setup. In comparing the left and the right panels of Fig. 4, we see 
that the addition of different nuclear targets improves the reconstruction for all benchmarks. 

In order to illustrate further the importance of accounting for astrophysical uncertainties 
for each of the models we consider, in Fig. 5 we show the bias in the reconstruction arising 
from assuming the wrong value for the outer slope c of the dark matter density profile 
and the scale radius r s when determining the velocity distribution through the Eddington 
formula. We have chosen the outer slope and the scale radius because those are two of the 
most significant parameters that impact the velocity distribution in our formalism [16, 34]. 
Though both halo models provide a good description to the galactic halo [16], it is clear 
that they imply very different underlying values for the WIMP properties. As above, the 
reconstruction is seen to be more biased in this case for lower mass WIMPs. We do find, 
however, that marginalizing over the outer slope and scale radius does reduce this bias, as 
shown by the upper sets of green (light) contours in Fig. 5. This simple procedure reduces 
the systematic bias by up to a factor of 5, thus significantly improving the accuracy of the 
reconstruction. 

In the context of our isotropic velocity distribution model, it is also possible for us to ask 
how well we can reconstruct the velocity distribution with an observed event rate distribution. 
In Fig. 6 we show the reconstructed velocity distribution assuming a 50 GeV mass WIMP 
and a ton-scale xenon detector. Overall we find that the reconstruction is unbiased at all 
velocities, so it is possible to determine the velocity distribution from direct detection data 
for these parameters, at least in parametric form. We also found that the marginal posterior 
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Figure 6. Reconstruction of the velocity distribution function for a ton-scale xenon detector and 
a 50 GeV WIMP. The inner and outer bands indicate the 1 and 2a ranges around the mean of the 
marginal posterior, respectively. The true velocity distribution is shown by the solid line. 

distribution at each velocity value is fairly close to Gaussian. Note that the measured event 
rate does not directly probe the velocity distribution below the minimum velocity to scatter 
at a given energy in xenon. Therefore, the parametric reconstruction presented here relies 
somewhat on the specific class of models we assumed. However, it is encouraging to note 
that, in this context, the reconstruction of the velocity distribution is unbiased. 

6 Conclusion 

In this paper we have addressed the question of the bias in the determination of the WIMP 
mass and spin-independent cross-section for plausible models of the velocity distribution 
function. For the two specific models that we consider, the standard Maxwell-Boltzmann 
distribution and isotropic, spherically symmetric model determined from the Eddington for- 
mula, we have shown that assuming the wrong velocity distribution model only leads to 
moderate bias in the determination of WIMP parameters. This can be further mitigated by 
using astrophysical data to constrain galactic model parameters and properly marginalizing 
over them. In addition to the different parameterizations used for the velocity distribution 
function in our modeling, we use kinematic data to marginalize over halo model parameters in 
determining the velocity distribution. Thus this paper represents an important step towards 
connecting WIMP event rates to a full galactic halo model in a theoretically self-consistent 
framework. 

While the focus of our analysis has been on two particular models of the velocity distri- 
bution, our method can be applied to analytic models that depend on both the energy and 
the angular momentum [20]. Though these models add degrees of freedom to the galactic 
model analysis, it is not clear whether they are a better description of our dark matter halo. 
However, anisotropic models do admit interesting theoretical phenomenology, allowing for 
potential model-independent extraction of the anisotropy of the WIMP velocity distribution 
[49]. 
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This paper focused on the simplest elastic scattering dark matter models for detectors 
without directional sensitivity. In the future it will also be interesting to consider directional 
signals. Initial analysis along these lines have been undertaken assuming phenomenological 
models of the velocity distribution [50]. Connecting these distributions to both galactic halo 
model parameters and to the dark matter distribution in the Galaxy, as we have done here 
for isotropic models, will likely provide more rigorous constraints on the three-dimensional 
WIMP distribution [51]. Further, it will be interesting to apply our methods to annual 
modulation signals [52]. Also, in the future it will be important to apply the method that 
we have outlined to more general parameterizations of the elastic scattering cross-section 
and spin-dependent cross-sections [53, 54]. In fact, spin-dependent couplings induce different 
dependencies of the event rate on the nuclear recoil energy and can in many cases lead to 
stronger bounds than spin-independent scattering [55]. 

A general conclusion of our study is that the use of the standard Maxwellian velocity 
distribution leads to a robust determination of the WIMP mass from direct detection exper- 
iments, as long as astrophysical uncertainties are kept into account through marginalization 
over galactic model parameters. 
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